#pragma once

//
// alt ("algebraic template library")
//
//written by Ryo Hashimoto(Apr 17, 2008)

namespace alt{

    typedef std::complex<double> complex_t;

    template<typename ValueType,size_t N>
    class algebraic_equation :
        boost::operators<algebraic_equation<ValueType,N>,ValueType >
	{
    public:
        typedef ValueType value_type;
        static const size_t order = N-1;
        typedef std::vector<value_type> solution_container;

        typedef algebraic_equation<value_type,N> this_type;
    private:
        value_type m_a[N];
        
        //operator implementation detail
        template<bool this_type_is_result_type,size_t N_r>
        struct add_impl{
            typedef this_type lhs_type;
            typedef algebraic_equation<value_type,N_r> rhs_type;
            typedef rhs_type result_type;

            result_type operator()(const lhs_type& lhs,const rhs_type& rhs)const
            {
                result_type res(rhs);
                for(size_t i=0;i!=N;++i){
                    res[i] += lhs[i];
                }
                return res;
            }
        };
        template<size_t N_r>
        struct add_impl<true,N_r>{
            typedef this_type lhs_type;
            typedef algebraic_equation<value_type,N_r> rhs_type;
            typedef this_type result_type;

            result_type operator()(const lhs_type& lhs,const rhs_type& rhs)const
            {
                result_type res(lhs);
                for(size_t i=0;i!=N_r;++i){
                    res[i] += rhs[i];
                }
                return res;
            }
        };

        template<bool this_type_is_result_type,size_t N_r>
        struct sub_impl{
            typedef this_type lhs_type;
            typedef algebraic_equation<value_type,N_r> rhs_type;
            typedef rhs_type result_type;
            
            result_type operator()(const lhs_type& lhs,const rhs_type& rhs)const
            {
                result_type res;
                for(size_t i=0;i!=N_r;++i){
                    res[i] = -rhs[i];
                }
                for(size_t i=0;i!=N;++i){
                    res[i] += lhs[i];
                }
                return res;
            }
        };
        template<size_t N_r>
        struct sub_impl<true,N_r>{
            typedef this_type result_type;
            typedef this_type lhs_type;
            typedef algebraic_equation<value_type,N_r> rhs_type;

            result_type operator()(const lhs_type& lhs,const rhs_type& rhs)const
            {
                result_type res(lhs);
                for(size_t i=0;i!=N_r;++i){
                    res[i] -= rhs[i];
                }
                return res;
            }
        };

    public:
        //constructor
        algebraic_equation()
        {
        }
        template<typename IterType>
        algebraic_equation(IterType begin,IterType end)
        {
            value_type* p = m_a;
            for(;begin!=end;++begin,++p){
                *p = *begin;
            }
        }
        //operators
        value_type& operator[](size_t index)
        {
            assert(index<N);
            return m_a[index];
        }
        const value_type& operator[](size_t index)const
        {
            assert(index<N);
            return m_a[index];
        }

        //add
        template<size_t M>
        algebraic_equation<value_type,((N>M)?N:M)>
            operator+(const algebraic_equation<value_type,M> &r)const
        {
            const size_t N_result = (N>M)?N:M;
            typedef algebraic_equation<value_type,N_result> result_type;
            return add_impl<N==N_result,M>()(*this,r);
        }
        this_type& operator+=(const value_type& r)
        {
            (*this)[0] += r;
            return *this;
        }

        //subtract
        template<size_t M>
        algebraic_equation<value_type,((N>M)?N:M)>
            operator-(const algebraic_equation<value_type,M> &r)const
        {
            const size_t N_result = (N>M)?N:M;
            typedef algebraic_equation<value_type,N_result> result_type;
            return sub_impl<N==N_result,M>()(*this,r);
        }
        this_type& operator-=(const value_type& r)
        {
            (*this)[0] -= r;
            return *this;
        }

        //multiply
        template<size_t M>
        algebraic_equation<value_type,N+M-1>
            operator*(const algebraic_equation<value_type,M> &r)const
        {
            const size_t N_result = N+M-1;
            typedef algebraic_equation<value_type,N_result> result_type;

            result_type res;

            for(size_t i=0;i!=N_result;++i)res[i]=0;

            for(size_t i=0;i!=N;++i){
                for(size_t j=0;j!=M;++j){
                    res[i+j] += (*this)[i]*r[j];
                }
            }

            return res;
        }
        this_type& operator*=(const value_type& r)
        {
            for(size_t i=0;i!=N;++i){
                (*this)[i] *= r;
            }
            return *this;
        }

        //divide
        this_type& operator/=(const value_type& r)
        {
            for(size_t i=0;i!=N;++i){
                (*this)[i] /= r;
            }
            return *this;
        }

        //equal
        template<size_t M>
        algebraic_equation<value_type,((N>M)?N:M)>
            operator==(const algebraic_equation<value_type,M> &r)const
        {
            const size_t N_result = (N>M)?N:M;
            typedef algebraic_equation<value_type,N_result> result_type;

            return (*this)-r;
        }
        this_type operator==(const value_type& r)
        {
            return *this-r;
        }

        //output
        friend std::ostream& operator<<(std::ostream& s,const this_type& eq)
        {
            for(size_t i=0;i!=N;++i){
                s << eq[N-1-i] << ((i<(N-1))?",":"");
            }
            return s;
        }
    };

    
    //power
    template<size_t M,size_t N,typename ValueType>
    struct pow_impl{
        algebraic_equation<ValueType,(N-1)*M+1>
            operator()(const algebraic_equation<ValueType,N>& eq)const
        {
            return eq*pow_impl<M-1,N,ValueType>()(eq);
        }
    };
    template<size_t N,typename ValueType>
    struct pow_impl<0,N,ValueType>{
        algebraic_equation<ValueType,1> 
            operator()(const algebraic_equation<ValueType,N>& eq)
        {
            typedef algebraic_equation<ValueType,1> result_type;
            result_type res;
            res[0] = 1;
            return res;
        }
    };
    template<size_t M,size_t N,typename ValueType>
    algebraic_equation<ValueType,(N-1)*M+1>
        pow(const algebraic_equation<ValueType,N>& eq)
    {
        return pow_impl<M,N,ValueType>()(eq);
    }
	
    //utility
	template<typename ValueType>
	class variant_impl : public algebraic_equation<ValueType,2>
	{
	public:
		typedef algebraic_equation<ValueType,2> super;
		variant_impl()
		{
			(*this)[1] = 1;
			(*this)[0] = 0;
		}
	};
	typedef variant_impl<complex_t> variant;

    //get solution
    template<size_t N,typename ValueType>//declaration
	typename variant_impl<ValueType>::solution_container get_solutions(const algebraic_equation<ValueType,N>& eq)
    {
    }

    template<typename ValueType>//linear equation
    typename variant_impl<ValueType>::solution_container
        get_solutions(const algebraic_equation<ValueType,2>& eq)
    {
		typedef ValueType value_type;
        typedef algebraic_equation<value_type,2> EqType;

        typename variant_impl<value_type>::solution_container res(EqType::order);

        res[0] = -eq[0]/eq[1];

        return res;
    }

    template<typename ValueType>//quadratic equation
    typename variant_impl<ValueType>::solution_container
        get_solutions(const algebraic_equation<ValueType,3>& eq)
    {
		typedef ValueType value_type;
        typedef algebraic_equation<value_type,3> EqType;

        const value_type& a = eq[2];
        const value_type& b = eq[1];
        const value_type& c = eq[0];

        const ValueType D = sqrt(b*b - value_type(4)*a*c);

        typename variant_impl<value_type>::solution_container res(EqType::order);

        res[0] = (-b+D)/(value_type(2)*a);
        res[1] = (-b-D)/(value_type(2)*a);

        return res;
    }

    //cubic equation
    variant_impl<complex_t>::solution_container
        get_solutions(const algebraic_equation<complex_t,4>& eq)
    {
        typedef complex_t value_type;
        typedef algebraic_equation<value_type,4> EqType;

        const value_type& a = eq[3];

        const value_type A2 = eq[2]/a;
        const value_type A1 = eq[1]/a;
        const value_type A0 = eq[0]/a;

        const value_type omega = value_type(-1.0/2,sqrt(3.0)/2);

        const value_type b = -(value_type(-2)*A2*A2*A2 + value_type(9)*A1*A2 - value_type(27)*A0);
        const value_type _c = (A2*A2 - value_type(3)*A1);
        const value_type c = _c*_c*_c;

        const value_type D = sqrt(b*b - value_type(4)*a*c);

        const value_type s1_3 = (-b+D)/value_type(2);
        const value_type s2_3 = (-b-D)/value_type(2);

        const value_type s0 = -A2;
        const value_type s1 = pow(s1_3,1.0/3);
        const value_type s2 = pow(s2_3,1.0/3);

        variant_impl<value_type>::solution_container res(EqType::order);

        res[0] = (s0 + s1 + s2)/value_type(3);
        res[1] = (s0 + omega*omega*s1 + omega*s2)/value_type(3);
        res[2] = (s0 + omega*s1 + omega*omega*s2)/value_type(3);

        return res;
    }
}
